{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "approximate-imperial",
   "metadata": {},
   "source": [
    "# 简单理解最初步的 ONIOM 能量计算"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "according-width",
   "metadata": {},
   "source": [
    "> 创建时间：2021-06-10"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sufficient-shanghai",
   "metadata": {},
   "source": [
    "在这篇文档中，我们会使用 PySCF 实现 ONIOM 计算的最初步的结果，即单点能计算。ONIOM 方法的原始文献之一是 Dapprich, Frisch et al. [^Dapprich-Frisch.JMST.1999]。\n",
    "\n",
    "由于**单点能本身是没有实际物理价值的**，因此相对能量、或者光谱性质、几何构型才是真正有价值的输出量。但得到这些输出量会比较复杂，特别是大多数 (除了 Gaussian 外) 以代价较大的第一性为卖点的量化程序，对半经验方法的光谱信息的支持并不好。我们并不对这些问题作更细致的讨论。正因为 Gaussian 对 ONIOM 的各种光谱数据、结构性质和溶剂化的支持非常充分，而其它软件在没有特化程序时似乎只能计算单点能；因此这些特色使得 Gaussian 确实地成为适合药物设计的软件之一。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "stock-screen",
   "metadata": {},
   "source": [
    "尽管 PySCF 确实可以给出 ONIOM 的计算结果，但为了方便地得到 ONIOM 的分块信息，我们仍然还需要 Gaussian 的支持。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "satisfied-significance",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyscf import gto, scf, mp, semiempirical\n",
    "import numpy as np\n",
    "\n",
    "HARTREE2KCAL = semiempirical.mopac_param.HARTREE2KCAL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "living-egypt",
   "metadata": {},
   "source": [
    "## 简单问题：三氟代乙醛的 2-layer 计算"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "persistent-technical",
   "metadata": {},
   "source": [
    "该问题的 Gaussian 输入卡置于 {download}`2-layer.gjf` 中。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "narrative-stretch",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#p ONIOM(MP2(Full)/6-31G:HF/6-31G)=InputFiles NoSymm\r\n",
      " \r\n",
      "2-layer ONIOM, modified from Gaussian keyword list\r\n",
      " \r\n",
      "0 1 0 1 0 1\r\n",
      "  F     -1.041506214819     0.000000000000    -2.126109488809 L\r\n",
      "  F     -2.033681935634    -1.142892069126    -0.412218766901 L\r\n",
      "  F     -2.033681935634     1.142892069126    -0.412218766901 L\r\n",
      "  C     -1.299038105677     0.000000000000    -0.750000000000 L H 5 0.7 0.8\r\n",
      "  C      0.000000000000     0.000000000000     0.000000000000 H\r\n",
      "  H      0.000000000000     0.000000000000     1.100000000000 H\r\n",
      "  O      1.125833024920     0.000000000000    -0.650000000000 H\r\n",
      "\r\n"
     ]
    }
   ],
   "source": [
    "! cat 2-layer.gjf"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "challenging-edinburgh",
   "metadata": {},
   "source": [
    "计算完毕后，最关键的 Gaussian 的输出是\n",
    "\n",
    "```\n",
    "ONIOM: gridpoint  1 method:  low   system:  model energy:  -113.796286020376\n",
    "ONIOM: gridpoint  2 method:  high  system:  model energy:  -114.023466233349\n",
    "ONIOM: gridpoint  3 method:  low   system:  real  energy:  -449.289045248409\n",
    "ONIOM: extrapolated energy =    -449.516225461381\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fresh-composition",
   "metadata": {},
   "source": [
    "我们可以知道，最终能量的计算是通过下式完成的：\n",
    "\n",
    "$$\n",
    "E_\\mathrm{ONIOM2} = E_\\mathrm{low} (\\mathrm{real}) + E_\\mathrm{high} (\\mathrm{model}) - E_\\mathrm{low} (\\mathrm{model})\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "taken-straight",
   "metadata": {},
   "source": [
    "对于 2-layer 的 ONIOM 计算，它要细分成三个计算任务：\n",
    "\n",
    "| 任务序号 | 计算级别 (method) | 计算体系 (system) | 能量 / Hartree |\n",
    "|--|--:|--:|--:|\n",
    "| 1 | HF/6-31G (low) | (H)CHO (model) | -113.796286020376 |\n",
    "| 2 | MP2(Full)/6-31G (high) | (H)CHO (model) | -114.023466233349 |\n",
    "| 3 | HF/6-31G (low) | CF3-CHO (real) | -449.289045248409 |"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "severe-holocaust",
   "metadata": {},
   "source": [
    "注意到我们指定的模型层 (model) 是最后三个原子；但实际上，由于模型层 (model) 与全局层 (real) 之间有键相互作用 (第 4 个 real 层甲基碳原子与第 5 个 model 层醛碳原子)，因此不能简单粗暴地在此处断键。\n",
    "\n",
    "上面的 Gaussian 输入卡会在计算模型层时引入氢原子；其引入并非直接将第 4 个碳原子替换为氢，而同时还要缩放键长度。在上述表格中，人为补上去的氢原子也由括号作标识。`H 5 0.7 0.8` 即指在低计算级别 (low) 缩放到 0.7 倍，而在高计算级别 (high) 缩放到 0.8 倍。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "drawn-graduation",
   "metadata": {},
   "source": [
    "### 任务 1：$E_\\mathrm{low} (\\mathrm{model})$ 计算"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "wireless-phone",
   "metadata": {},
   "source": [
    "低计算级别 (low) 与高计算级别 (high) 都需要计算模型层 (model) 的分子片段。但由于氢原子添补方式的不同，使得我们需要分别定义这两个分子，及其后续计算。我们首先给出 high 级别的计算结果。\n",
    "\n",
    "Low 级别下，引入的氢原子键长是对应 4 号碳与 5 号碳的 0.7 倍。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "foreign-brisbane",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<pyscf.gto.mole.Mole at 0x7f277b9a56a0>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mol_model_low = gto.Mole()\n",
    "mol_model_low.atom = \"\"\"\n",
    "  H     -0.90932667         0.                -0.525     \n",
    "  C      0.000000000000     0.000000000000     0.000000000000 \n",
    "  H      0.000000000000     0.000000000000     1.100000000000 \n",
    "  O      1.125833024920     0.000000000000    -0.650000000000 \n",
    "\"\"\"\n",
    "mol_model_low.basis = \"6-31G\"\n",
    "mol_model_low.verbose = 0\n",
    "mol_model_low.build()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "moderate-manhattan",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-113.79628602351522"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mf_model_low = scf.RHF(mol_model_low).run()\n",
    "eng_model_low = mf_model_low.e_tot\n",
    "eng_model_low"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "usual-narrow",
   "metadata": {},
   "source": [
    "### 任务 2：$E_\\mathrm{high} (\\mathrm{model})$ 计算"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "affiliated-burden",
   "metadata": {},
   "source": [
    "High 级别下，引入的氢原子键长是对应 4 号碳与 5 号碳的 0.8 倍。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "harmful-symphony",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<pyscf.gto.mole.Mole at 0x7f276ae590a0>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mol_model_high = gto.Mole()\n",
    "mol_model_high.atom = \"\"\"\n",
    "  H     -1.03923048         0.                -0.6\n",
    "  C      0.000000000000     0.000000000000     0.000000000000 \n",
    "  H      0.000000000000     0.000000000000     1.100000000000 \n",
    "  O      1.125833024920     0.000000000000    -0.650000000000 \n",
    "\"\"\"\n",
    "mol_model_high.basis = \"6-31G\"\n",
    "mol_model_high.verbose = 0\n",
    "mol_model_high.build()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "advanced-functionality",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-114.02346625303319"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mf_model_high = mp.MP2(mol_model_high).run()\n",
    "eng_model_high = mf_model_high.e_tot\n",
    "eng_model_high"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "spiritual-webcam",
   "metadata": {},
   "source": [
    "### 任务 3：$E_\\mathrm{low} (\\mathrm{real})$ 计算"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "comparable-collar",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<pyscf.gto.mole.Mole at 0x7f276ae597c0>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mol_real = gto.Mole()\n",
    "mol_real.atom = \"\"\"\n",
    "  F     -1.041506214819     0.000000000000    -2.126109488809 \n",
    "  F     -2.033681935634    -1.142892069126    -0.412218766901 \n",
    "  F     -2.033681935634     1.142892069126    -0.412218766901 \n",
    "  C     -1.299038105677     0.000000000000    -0.750000000000 \n",
    "  C      0.000000000000     0.000000000000     0.000000000000 \n",
    "  H      0.000000000000     0.000000000000     1.100000000000 \n",
    "  O      1.125833024920     0.000000000000    -0.650000000000 \n",
    "\"\"\"\n",
    "mol_real.basis = \"6-31G\"\n",
    "mol_real.verbose = 0\n",
    "mol_real.build()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "optional-discretion",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-449.28904521819294"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mf_real_low = scf.RHF(mol_real).run()\n",
    "eng_real_low = mf_real_low.e_tot\n",
    "eng_real_low"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "beneficial-going",
   "metadata": {},
   "source": [
    "### 能量的统合"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "delayed-reproduction",
   "metadata": {},
   "source": [
    "回顾 2-layer ONIOM 的能量统合方式：\n",
    "\n",
    "$$\n",
    "E_\\mathrm{ONIOM2} = E_\\mathrm{low} (\\mathrm{real}) + E_\\mathrm{high} (\\mathrm{model}) - E_\\mathrm{low} (\\mathrm{model})\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "arbitrary-affairs",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-449.5162254477109"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eng_real_low + eng_model_high - eng_model_low"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "improving-entrepreneur",
   "metadata": {},
   "source": [
    "Gaussian 的结果是 -449.516225461381。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "finished-emphasis",
   "metadata": {},
   "source": [
    "## 较复杂问题：丙醛的 3-layer 计算"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "whole-phenomenon",
   "metadata": {},
   "source": [
    "该问题的 Gaussian 输入卡置于 {download}`3-layer.gjf` 中。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "proved-monday",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#p ONIOM(MP2(Full)/6-311g*:HF/6-31g:MINDO3)=InputFiles NoSymm\r\n",
      "\r\n",
      "3-layer ONIOM, modified from Gaussian test job 0679\r\n",
      "\r\n",
      "   0   1   0   1   0   1   0   1   0   1   0   1   0   1\r\n",
      " C  -0.006049274275    0.000000000000    0.066754956170 H\r\n",
      " O   0.011403425950    0.000000000000    1.308239478983 H\r\n",
      " H   0.944762558657    0.000000000000   -0.507359536461 H\r\n",
      " C  -1.307562483867    0.000000000000   -0.766510748030 M H 1 0.723886 0.723886 0.723886\r\n",
      " C  -1.047480751885    0.000000000000   -2.301387120377 L H 4 0.723886 0.723886 0.723886\r\n",
      " H  -1.903669606697   -0.885256630266   -0.468844831106 M\r\n",
      " H  -1.903669606697    0.885256630266   -0.468844831106 M\r\n",
      " H  -1.988817319373    0.000000000000   -2.842389774687 L\r\n",
      " H  -0.482972255230    0.881286097766   -2.591806824941 L\r\n",
      " H  -0.482972255230   -0.881286097766   -2.591806824941 L\r\n",
      "\r\n",
      "\r\n"
     ]
    }
   ],
   "source": [
    "! cat 3-layer.gjf"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "quality-bonus",
   "metadata": {},
   "source": [
    "计算完毕后，最关键的 Gaussian 的输出是\n",
    "\n",
    "```\n",
    "ONIOM: gridpoint  1 method:  low   system:  model energy:    -0.033309689782\n",
    "ONIOM: gridpoint  2 method:  med   system:  model energy:  -113.805007046900\n",
    "ONIOM: gridpoint  3 method:  low   system:  mid   energy:    -0.059535553265\n",
    "ONIOM: gridpoint  4 method:  high  system:  model energy:  -114.255323473041\n",
    "ONIOM: gridpoint  5 method:  med   system:  mid   energy:  -152.836036428501\n",
    "ONIOM: gridpoint  6 method:  low   system:  real  energy:    -0.068015765875\n",
    "ONIOM: extrapolated energy =    -153.294833067253\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "lasting-romantic",
   "metadata": {},
   "source": [
    "由于外推表达式比较复杂，原始文献也使用了简化的表达式作能量结果的表示 (需要注意，$E_1 = E_\\mathrm{low} (\\mathrm{model})$ 没有在能量统合公式中)：\n",
    "\n",
    "$$\n",
    "E_\\mathrm{ONIOM3} = E_6 - E_3 + E_5 - E_2 + E_4\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "unavailable-calibration",
   "metadata": {},
   "source": [
    "对于 3-layer 的 ONIOM 计算，它要细分成 5 个计算任务：\n",
    "\n",
    "| 任务序号 | 计算级别 (method) | 计算体系 (system) | 能量 / Hartree |\n",
    "|--|--:|--:|--:|\n",
    "| 2 | HF/6-31g (med) | (H)CHO (model) | -113.805007046900 |\n",
    "| 3 | MINDO/3 (low) | (H)CH2-CHO (mid) | -0.059535553265 |\n",
    "| 4 | MP2(Full)/6-311g* (high) | (H)CHO (model) | -114.255323473041 |\n",
    "| 5 | HF/6-31g (med) | (H)CH2-CHO (mid) | -152.836036428501 |\n",
    "| 6 | MINDO/3 (low) | CH3-CH2-CHO (real) | -0.068015765875 |"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "treated-variety",
   "metadata": {},
   "source": [
    "### 使用 Gaussian 查看每个计算任务的分子构型"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "proud-unemployment",
   "metadata": {},
   "source": [
    "一旦体系增大，模型 (model)、中间 (mid/intermediate) 与全局 (real) 层之间的断键数目增多，这时人为地确定具体的分子构型就会变得困难了。\n",
    "\n",
    "一种比较简单粗暴的方式是使用 Gaussian ONIOM 关键词的 `OnlyInputFiles` 选项；这样就可以在不执行 ONIOM 具体计算的情况下，把每个需要计算的分子片段信息打印出来。如果加入关键词 `InputFiles`，那么可以同时打印分子片段以及计算 ONIOM。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "celtic-vermont",
   "metadata": {},
   "source": [
    "我们以第 5 个格点为例 (中层 (mid/intermediate)，中等级算量 (med))。其在 Gaussian 中的输出是 (略去 Gaussian 推断的成键关系)\n",
    "\n",
    "```\n",
    "ONIOM: generating point  5 -- med level on mid system.\n",
    "\n",
    "--------------------------------------------------------------------------------\n",
    "#P Test IOp(2/15=1,5/32=2,5/38=1) HF/6-31G\n",
    "\n",
    "3-layer ONIOM, modified from Gaussian test job 0679\n",
    "Point  5 -- med level on mid system.\n",
    "\n",
    "0,1\n",
    " C                                                -0.006049274275      0.000000000000      0.066754956170\n",
    " O                                                 0.011403425950      0.000000000000      1.308239478983\n",
    " H                                                 0.944762558657      0.000000000000     -0.507359536461\n",
    " C                                                -1.307562483867      0.000000000000     -0.766510748030\n",
    " H(Iso=12)                                        -1.119292959229      0.000000000000     -1.877586265703\n",
    " H                                                -1.903669606697     -0.885256630266     -0.468844831106\n",
    " H                                                -1.903669606697      0.885256630266     -0.468844831106\n",
    " Bq-#1(Iso=1.00782504,Spin=1,ZNuc=1.,GFac=2.792846)         -1.988817319373      0.000000000000     -2.842389774687\n",
    " Bq-#1(Iso=1.00782504,Spin=1,ZNuc=1.,GFac=2.792846)         -0.482972255230      0.881286097766     -2.591806824941\n",
    " Bq-#1(Iso=1.00782504,Spin=1,ZNuc=1.,GFac=2.792846)         -0.482972255230     -0.881286097766     -2.591806824941\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "disciplinary-revision",
   "metadata": {},
   "source": [
    "我们可以看出第 5 号位的氢原子实际上是替代了低层 (low) 的碳原子。之所以要设置为氢原子 12 质量的同位素，是为了跟进一步的分子力与频率分析作准备；在单纯讨论能量时不需要考虑原子质量的问题。同时，最后的三个低层 (low) 原子都被设置为 `Bq` 原子，即虚原子。我们实际需要带入能量计算的原子就是前 7 个原子了。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "electoral-chorus",
   "metadata": {},
   "source": [
    "### 任务 2：$E_2 = E_\\mathrm{med} (\\mathrm{model})$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "cellular-revision",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<pyscf.gto.mole.Mole at 0x7f276ae59a30>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mol_2 = gto.Mole()\n",
    "mol_2.atom = \"\"\"\n",
    "C    -0.006049274275      0.000000000000      0.066754956170\n",
    "O     0.011403425950      0.000000000000      1.308239478983\n",
    "H     0.944762558657      0.000000000000     -0.507359536461\n",
    "H    -0.948196465514      0.000000000000     -0.536434421381\n",
    "\"\"\"\n",
    "mol_2.basis = \"6-31G\"\n",
    "mol_2.verbose = 0\n",
    "mol_2.build()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "middle-boston",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-113.80500704910236"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mf_2 = scf.RHF(mol_2).run()\n",
    "eng_2 = mf_2.e_tot\n",
    "eng_2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "lonely-offer",
   "metadata": {},
   "source": [
    "### 任务 3：$E_3 = E_\\mathrm{low} (\\mathrm{mid})$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "oriented-grave",
   "metadata": {},
   "source": [
    "需要注意，Gaussian 在半经验方法中输出的能量并非接近于单点能，而是接近于原子化能。因此，在使用 PySCF 时尽量不要直接用半经验的 `e_tot` 变量作结果输出。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "atmospheric-guidance",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<pyscf.gto.mole.Mole at 0x7f276ae59a60>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mol_3 = gto.Mole()\n",
    "mol_3.atom = \"\"\"\n",
    "C     -0.006049274275      0.000000000000      0.066754956170\n",
    "O      0.011403425950      0.000000000000      1.308239478983\n",
    "H      0.944762558657      0.000000000000     -0.507359536461\n",
    "C     -1.307562483867      0.000000000000     -0.766510748030\n",
    "H     -1.119292959229      0.000000000000     -1.877586265703\n",
    "H     -1.903669606697     -0.885256630266     -0.468844831106\n",
    "H     -1.903669606697      0.885256630266     -0.468844831106\n",
    "\"\"\"\n",
    "mol_3.verbose = 0\n",
    "mol_3.build()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "celtic-sullivan",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.059543662948324035"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mf_3 = semiempirical.MINDO3(mol_3).run()\n",
    "eng_3 = mf_3.e_heat_formation / HARTREE2KCAL\n",
    "eng_3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "operating-consultancy",
   "metadata": {},
   "source": [
    "这与 Gaussian 的结果有略微差别，但差别在 5e-3 kcal/mol 上，我们可以认为这时可以忽略的差距了。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "through-controversy",
   "metadata": {},
   "source": [
    "### 任务 4：$E_4 = E_\\mathrm{high} (\\mathrm{model})$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "embedded-manitoba",
   "metadata": {},
   "source": [
    "需要注意，即使这个分子与 $E_2 = E_\\mathrm{med} (\\mathrm{model})$ 所使用的分子相同 (因为使用了相同的断键补氢系数 $g = 0.723886$)；但 $E_2$ 的基组是中等级 (med) 的 6-31G，而 $E_4$ 则是高等级 (high) 的 6-311G*。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "convertible-consideration",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<pyscf.gto.mole.Mole at 0x7f276ae598e0>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mol_4 = gto.Mole()\n",
    "mol_4.atom = \"\"\"\n",
    "C    -0.006049274275      0.000000000000      0.066754956170\n",
    "O     0.011403425950      0.000000000000      1.308239478983\n",
    "H     0.944762558657      0.000000000000     -0.507359536461\n",
    "H    -0.948196465514      0.000000000000     -0.536434421381\n",
    "\"\"\"\n",
    "mol_4.basis = \"6-311G*\"\n",
    "mol_4.verbose = 0\n",
    "mol_4.build()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "liberal-cargo",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-114.25532346425724"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mf_4 = mp.MP2(mol_4).run()\n",
    "eng_4 = mf_4.e_tot\n",
    "eng_4"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "considered-belief",
   "metadata": {},
   "source": [
    "### 任务 5：$E_5 = E_\\mathrm{med} (\\mathrm{mid})$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "previous-logistics",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<pyscf.gto.mole.Mole at 0x7f276ae77670>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mol_5 = gto.Mole()\n",
    "mol_5.atom = \"\"\"\n",
    "C     -0.006049274275      0.000000000000      0.066754956170\n",
    "O      0.011403425950      0.000000000000      1.308239478983\n",
    "H      0.944762558657      0.000000000000     -0.507359536461\n",
    "C     -1.307562483867      0.000000000000     -0.766510748030\n",
    "H     -1.119292959229      0.000000000000     -1.877586265703\n",
    "H     -1.903669606697     -0.885256630266     -0.468844831106\n",
    "H     -1.903669606697      0.885256630266     -0.468844831106\n",
    "\"\"\"\n",
    "mol_5.basis = \"6-31G\"\n",
    "mol_5.verbose = 0\n",
    "mol_5.build()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "bibliographic-fluid",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-152.83603642639747"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mf_5 = scf.RHF(mol_5).run()\n",
    "eng_5 = mf_5.e_tot\n",
    "eng_5"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "complex-promise",
   "metadata": {},
   "source": [
    "### 任务 6：$E_6 = E_\\mathrm{low} (\\mathrm{real})$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "female-momentum",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<pyscf.gto.mole.Mole at 0x7f276ae77d90>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mol_6 = gto.Mole()\n",
    "mol_6.atom = \"\"\"\n",
    "C     -0.006049274275      0.000000000000      0.066754956170\n",
    "O      0.011403425950      0.000000000000      1.308239478983\n",
    "H      0.944762558657      0.000000000000     -0.507359536461\n",
    "C     -1.307562483867      0.000000000000     -0.766510748030\n",
    "C     -1.047480751885      0.000000000000     -2.301387120377\n",
    "H     -1.903669606697     -0.885256630266     -0.468844831106\n",
    "H     -1.903669606697      0.885256630266     -0.468844831106\n",
    "H     -1.988817319373      0.000000000000     -2.842389774687\n",
    "H     -0.482972255230      0.881286097766     -2.591806824941\n",
    "H     -0.482972255230     -0.881286097766     -2.591806824941\n",
    "\"\"\"\n",
    "mol_6.verbose = 0\n",
    "mol_6.build()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "subject-remains",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.06801695739152386"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mf_6 = semiempirical.MINDO3(mol_6).run()\n",
    "eng_6 = mf_6.e_heat_formation / HARTREE2KCAL\n",
    "eng_6"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "drawn-chain",
   "metadata": {},
   "source": [
    "### 能量的统合"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "brilliant-highlight",
   "metadata": {},
   "source": [
    "回顾 3-layer ONIOM 的能量统合方式：\n",
    "\n",
    "$$\n",
    "E_\\mathrm{ONIOM3} = E_6 - E_3 + E_5 - E_2 + E_4\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "third-variable",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-153.29482613599555"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eng_6 - eng_3 + eng_5 - eng_2 + eng_4"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "minute-balloon",
   "metadata": {},
   "source": [
    "Gaussian 的结果是 -153.294833067253 Hartree。我们使用 PySCF 给出的计算结果与 Gaussian 相差 4e-3 kcal/mol。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "opposite-cotton",
   "metadata": {},
   "source": [
    "[^Dapprich-Frisch.JMST.1999]: Dapprich, S.; Komáromi, I.; Byun, K. S.; Morokuma, K.; Frisch, M. J. A New ONIOM Implementation in Gaussian98. Part I. the Calculation of Energies, Gradients, Vibrational Frequencies and Electric Field Derivatives. *J. Mol. Struct. THEOCHEM* **1999**, *461-462*, 1–21. doi: [10.1016/s0166-1280(98)00475-8](https://doi.org/10.1016/s0166-1280(98)00475-8)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "194.2px"
   },
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
